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ABSTRACT 

We have searched for planetary companions around 800 pulsars monitored at the Jodrell Bank Observatory, with both circular 
and eccentric orbits of periods between 20 days and 17 years and inclination-dependent planetary masses from 1074 to 100 Mọ. 
Using a Bayesian framework, we simultaneously model pulsar timing parameters and a stationary noise process with a power-law 
power spectral density. We put limits on the projected masses of any planetary companions, which reach as low as 1/100th of 
the mass of the Moon (~ 1074 Me). We find that two-thirds of our pulsars are highly unlikely to host any companions above 
2 — 8 Moe. Our results imply that fewer than 0.5% of pulsars could host terrestrial planets as large as those known to orbit PSR 
B1257+12 (~ 4Mga); however, the smaller planet in this system (~ 0.02 Mẹ) would be undetectable in 95% of our sample, 
hidden by both instrumental and intrinsic noise processes, although it is not clear if such tiny planets could exist in isolation. We 
detect significant periodicities in 15 pulsars, however we find that intrinsic quasi-periodic magnetospheric effects can mimic the 
influence of a planet, and for the majority of these cases we believe this to be the origin of the detected periodicity. Notably, we 
find that the highly periodic oscillations in PSR B0144+59 are correlated with changes in the pulse profile and therefore can be 
attributed to magnetospheric effects. We believe the most plausible candidate for planetary companions in our sample is PSR 
J20074+3120. 
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detection; methods: data analysis 


1 INTRODUCTION 


The rapid and regular pulsations observed from pulsars allows us to 
track their rotation to a small fraction of the spin period over long 
timescales. This technique of pulsar timing allows us to use them as 
tools for a wide range of physics and astrophysics applications, and 
makes our observations very sensitive to perturbations due to orbital 
motion, even from Earth or Moon mass planets (Thorsett & Phillips 
1992; Cordes 1993). 

Notably, the first extrasolar planets were discovered around the pul- 
sar PSR B1257+12 (Wolszczan & Frail 1992; Wolszczan 1994), and 
are still among the lowest mass extrasolar planets known with masses 
of 0.020(2) Mg, 4.3(2) Me and 3.9(2) Mg (Konacki & Wolszczan 
2003). Since then, almost 5000 exoplanets have been discovered by 
various methods (see e.g. The Extrasolar Planets Encyclopaedia!). 
However, in contrast with this considerable increase in the number of 
exoplanets discovered around other types of stars, only five other pul- 
sars have been confirmed to have planetary-mass companions. Of the 
pulsars with planetary-mass companions, four are ultra-low mass car- 
bon white dwarfs, or “diamond planets’, around pulsars J1719—1438 
(Bailes et al. 2011), JO636+5128 (Stovall et al. 2014), J1311—3430 
(Romani et al. 2012; Pletsch et al. 2012) and J2322—2650 (Spiewak 
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et al. 2018), respectively. Lastly, pulsar B1620—26 is part of a triple 
system, including a super-Jupiter mass planet; however, the location 
of this system in the globular cluster M4 suggests the planet was cap- 
tured, rather than formed around the pulsar (Sigurdsson et al. 2003), 
a process unlikely to happen for pulsars outside of the high density of 
such a cluster. PSR B1257+12 therefore remains as a unique example 
of a system of approximately Earth-mass planets orbiting a neutron 
star. 


There have been several previous systematic searches for low-mass 
planetary companions around pulsars, none of which have reported 
any detections. Kerr et al. (2015, K15 hereafter) searched for planets 
around 151 (young) pulsars, a sample two orders of magnitude larger 
than previous efforts by Thorsett & Phillips (1992). K15 concluded 
that no planets more massive than 0.4 Mg and with periods of less 
than one year are found around their sampled pulsars. Most recently, 
Behrens et al. (2020) used Pulsar Timing Array quality data to set 
mass constraints as low as a Moon mass (~ 0.01 Mẹ) for possible 
planetary masses around 45 millisecond pulsars (MSPs), for periods 
between 7 and 2000 days, also reporting no detectable planets. 


The apparent rarity of systems like that of PSR B1257+12 may well 
be a consequence of the extreme conditions in which pulsars form. 
There have been several mechanisms proposed for the formation 
of planets around pulsars (see e.g. Podsiadlowski 1993; Phillips & 
Thorsett 1994). Of these, three different types of scenarios can be 
distinguished. 
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In some of the proposed scenarios, the planets had formed around 
a normal star and were then captured by a neutron star through the 
collision of the two stars (Podsiadlowski et al. 1991) or survived 
the subsequent evolution of the initial system towards a neutron star 
system (Bailes et al. 1991). The resulting system would consist of a 
normal pulsar with planetary companions in eccentric orbits. These 
types of events would be very rare as they require fine tuning of 
environment conditions for the planets to survive or for the collision 
to take place outside of a globular cluster. 


In the majority of the proposed scenarios, planets form around 
pulsars from a disk of material around the neutron star. In the simplest 
case, some of the mass ejected from the supernova in which the 
pulsar was formed does not fully escape the gravitational pull of the 
star and falls back into a disk (Bailes et al. 1991; Lin et al. 1991); 
subsequently, a normal pulsar, surrounded by relatively small mass 
planets in circular orbits, is expected. In other models, a more massive 
circumstellar disk is a product of the neutron star distrupting its binary 
companion. If this companion was a white dwarf (Podsiadlowski et al. 
1991; Van den Heuvel 1992), we might expect to form a millisecond 
pulsar with planets in circular orbits. 


A third type of mechanisms proposes that the planets are the final 
stage of evolution of binary millisecond pulsars. Depending on the 
model, the initial pulsar companion is either partially evaporated 
(Bailes et al. 1991; Krolik 1991; Rasio et al. 1992) or fully disrupted 
(Stevens et al. 1992) by the influence of the neutron star, creating 
either a singular planet from the remains of the original companion 
(such as a “diamond planet”), or a debris disk around the neutron 
star, which would then collapse to form planets in circular orbits. 


Overall, there are a large number of proposed formation paths 
of planets around pulsars, and therefore large scale searches of 
planetary-mass companions and their orbital parameters are crucial 
to constraining and determining the feasibility of various models. 
However, although the time of arrival (ToA) measurements of indi- 
vidual observations can be very precise, the detectability of planets 
around pulsars is also limited by the presence of so-called “timing 
noise” that manifests as a long-term red noise process in the rotation 
of the pulsar. Models of timing noise typically assume a stochastic 
stationary red noise process (e.g. Shannon & Cordes 2010), but there 
is also strong evidence for quasi-periodic switching between spin- 
down states, which has been observed to be correlated to pulse shape 
changes in some pulsars (Hobbs et al. 2010; Lyne et al. 2010). This 
presents a further challenge in searching for planetary companions, 
as it can not only mask binary signatures, but also mimic them. A 
careful investigation is therefore necessary, as there have been previ- 
ous claims of planet detections around pulsars (e.g. PSR B0329+54; 
Shabanova 1995), later proved to be caused by intrinsic pulsar timing 
noise (Konacki et al. 1999). 


Likewise, previous attempts have been made to explain two peri- 
odicities seen in PSR B1828—11 using a planetary model, but Bailes 
et al. (1993) found that modelling this behaviour using only two 
simple sinusoids was not enough to fully account for the observed 
variation in the ToAs. The quasi-periodic behaviour in B1828—11 
is now known to be of intrinsic magnetospheric origin (Lyne et al. 
2010; Stairs et al. 2019). To better understand the population of plan- 
etary systems in pulsars, it is therefore important to both look at a 
large sample of pulsars that have been observed for long time-spans, 
and also to carefully model the pulsar timing noise to identify any 
purely periodic signals arising from planetary companions. The pul- 
sar timing archive at the Jodrell Bank Observatory (JBO) is an ideal 
dataset for this project, consisting of regular observations of almost 
800 pulsars with the 76-m Lovell Telescope. This represents nearly a 
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third of all known pulsars, with observations spanning up to 50 years 
for each pulsar. 

In this work we apply Bayesian techniques for modelling timing 
noise to the available JBO data, and conduct the largest systematic 
search for pulsar planetary companions to date, as well as set limits 
to any existing population of planetary-mass bodies around pulsars. 
This paper is structured as follows. In Section 2, the observational 
properties are summarised. Section 3 introduces the parameters of a 
binary (planetary) system. The method and relevant pulsar parame- 
ters used in this analysis are further described in Section 4. The main 
results are presented and discussed in Section 5. 


2 OBSERVATIONS 


The dataset used in this work is composed of observations of approx- 
imately 800 pulsars from the Jodrell Bank pulsar timing database. 
All pulsars have observations from the 76-m Lovell Telescope. Ob- 
servations since April 2009 are recorded using a digital filterbank (a 
clone of the DFB described in Manchester et al. 2013), and processed 
using PSRCHIVE to form time and frequency averaged pulse profiles 
with a 200—400 MHz bandwidth centred at 1520 MHz (depending 
on RFI excision). RFI is excised by manual removal of small outlier 
segments of the band (“channels”) and portions of the observations 
(“sub-integrations”). ToAs are generated using PsRCHIVE with an 
analytic standard template derived from the sum of all DFB obser- 
vations. Prior to 2009 the data used a range of receivers and were 
recorded with several generations of analogue systems, processed us- 
ing legacy software to directly generate ToAs. For 50 pulsars, daily 
observations with the 25-m Mark II Telescope at Jodrell Bank were 
used to fill intervals when the Lovell Telescope was undergoing up- 
grades and repairs. The recording and processing setup is identical to 
that used for the Lovell observations. For more details on the previous 
receivers, see Hobbs et al. (2004). 

For 16 pulsars these data are supplemented with observations made 
using the 13-m “42-ft” telescope at the Jodrell Bank Observatory. 
These consist of daily observations since 1990 (depending on the 
pulsar) with a 5 MHz bandwidth centred at 610 MHz. Since 2011 
these observations used the “cobra2” coherent dedispersion backend. 
Additionally the earliest ToAs for 29 pulsars come from observations 
made by the NASA Deep Space Network (Downs & Reichley 1983; 
Downs & Krause-Polstorff 1986). 

In total the dataset consists of 800 pulsars, of which 104 are 
millisecond pulsars. In total there are ~ 30 000 Lovell ToAs (with 1/3 
from the DFB backend) covering ~ 17 000 years of pulsar rotational 
history, implying an average cadence of 20 days. However, it should 
be noted that the cadence varies significantly from pulsar to pulsar 
depending on the science case for any individual source. 

Observations do not record absolute phase, so an initial tim- 
ing solution is used to predict/estimate absolute phase. With suf- 
ficient observational cadence, this can be done unambiguously and a 
pulse number assigned to each ToA. Phase ambiguities arise at large 
glitches (instantaneous changes in spin-down), and these are left as 
free parameters as part of the pulsar timing model. Further details of 
the dataset and the study of the pulsar timing noise and astrometric 
parameters will be included in future publications. 


3 MODELLING THE INFLUENCE OF PLANETS 


When a pulsar is part of a binary system (either with a star or a 
planet), it revolves around the centre of mass of the system, moving 
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with respect to the observer on Earth. This influences the arrival 
time of the signal to the observer, and is known as Rømer delay. 
In the timing model, the Rømer delay is expressed as (Blandford & 
Teukolsky 1976) 


[ (cos E(t)—e) sinw+sin E(t) V1 — e? cos w], 
(1) 


where ż is the ToA; i is the inclination angle of the plane of the orbit 
in the sky; e is the eccentricity; apsp is the semi-major axis of the 
pulsar around the centre of mass, given by 


1/3 


APSR sini 
Ars (t) Sa 
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Gm p? 
4r? (m + Mpsr)2 


(2) 


aPSR = 


where G is the gravitational constant, Mpspr is the mass of the pulsar, 
m is the mass of the companion, and Py is the orbital period; w is the 
argument of the periapsis, i.e. the angle of the periapsis with respect 
to the plane of the sky; E(t) is the eccentric anomaly at time t and 
is related to the true anomaly, A(t), i.e. the angle of an object on 
the elliptical orbit with respect to the periapsis by (see e.g. Karttunen 
et al. 2007) 


e + cos Ar(t) 


E(t) = ——————.. 
cos EG) 1 + ecos A(t) 


(3) 

The mean anomaly (M (t)), which is the fraction of the orbit that 
has been completed at time f since the orbiting body passed periapsis 
at time fo, expressed as an angle between 0 and 27x, can also be 
expressed in terms of the eccentric anomaly by 


M(t)= = (t- to) = E(t) -e sin E(t). (4) 


If to is known, the two expressions of M (t) in Eq. 4 can be used to find 
E(t), which can then be used in Eq. 1 to calculate the Romer delay 
at each ToA. However, if Pp is shorter than the time span of the data, 
there are multiple (degenerate) solutions for fg, as there are multiple 
periapsis crossings in the observed time. To avoid this degeneracy, 
we considered a phase parameter (¢), taking values between 0 and 1, 
such that 


a a (5) 
1 

where we chose a reference time trgp = 55000 MJD. The value of & 
can be used to find a single tọ solution (which is the periapsis crossing 
closest to trap), as follows. Eq. 5 gives AT(tpERp) straightforwardly 
from ¢. Then, we compute cos E(tRgF) using Eq. 3 at £ = tREF 
and the eccentricity value e. Lastly, equating the two expressions for 
M(trgp) as in Eq. 4, and given E(trpr) from the previous step, e 
and Py, we can determine fo. 


$ 


4 THE ANALYSIS 
4.1 Method overview 


We implemented a Bayesian method to analyse pulsar timing data, 
by adapting the existing pulsar timing software ENTERPRISE (Ellis 
et al. 2019). This method consists of marginalising over parameters 
describing the deterministic timing model, and simultaneously fitting 
for white noise, red noise and a planet orbit. This approach is well- 
suited for a systematic search of planets around pulsars, for placing 
limits on the mass of any orbiting celestial bodies, and therefore for 
inferring statistically significant properties of the population of these 
planets. 


In practice, a function containing the Keplerian orbit parameters 
was added to the ENTERPRISE likelihood function, allowing for si- 
multaneous fitting with the existing functionality; this is referred 
in this work as the “extended ENTERPRISE” and integrated into the 
RUN_ENTERPRISE pulsar analysis toolset (Keith et al. 2022). The 
method used consists, broadly, of the following steps applied to the 
timing data of one pulsar: 


(i) The orbital period prior is divided into equally-spaced bins in 
log-period, and the extended ENTERPRISE is run for each of these bins, 
fitting for white noise, red noise, and planet orbital parameters (see 
Section 4.3). 

(ii) For each period bin, 3-0 “detections” in the mass posteriors 
(samples) are flagged for further investigation. To make our search 
for periodicities more thorough, and to counter the smearing effects 
of large period ranges, the mass posteriors in each of these period 
bins are divided into five equally-spaced period slices (“sub-bins”), 
and the 3-o check performed again for each of these slices. Note 
that, while these sub-bins are not guaranteed to each have enough 
mass samples to give reliable mass limits, we manually verify that 
any detections seen are genuine, and not due to a small number of 
mass samples. We therefore only use the mass posteriors in these 
sub-bins to check for detections, not to estimate mass limits. 

(iii) To further study the population of planets around pulsars, 
we also estimate limits for the mass of planets detectable with our 
method and data, as the 95-percentile in the mass posteriors in each 
period bin. 


4.2 Timing analysis 
4.2.1 Pulsar parameters 


There are a number of deterministic parameters that are known to 
affect the timing of most pulsars. For each pulsar in our analysis we 
fit for the spin frequency (v) and two frequency derivatives (Y and ¥); 
astrometric parameters (i.e. position, proper motion); orbital param- 
eters of known binary companions; and the clock offsets (“jumps”) 
between different telescope backend systems, if relevant. We then 
analytically marginalise over these parameters, using a design ma- 
trix computed by TEMpo2 (Edwards et al. 2006). Glitch parameters 
as given in Basu et al. (2021) were also used. 


4.2.2 Timing noise 


Pulsar timing is also affected by noise, which we model with a 
white-noise component, with parameters EFAC (Hobbs et al. 2006) 
and EQUAD (Liu et al. 2012), and a power-law red-noise compo- 
nent. We characterise the red noise by the index y and the log- 
amplitude log;g(Ayeq), as used in the typical power-law Fourier- 
domain Gaussian-process model of timing noise (see Lentati et al. 
2014; Van Haasteren et al. 2011), 


AZ =y 
a yr, (6) 


1272 \lyr7! 
where f is the frequency in the Fourier domain, and P(f) is the 
one-sided power spectral density. 


P(f)= 


4.3 Planet fitting 


The set of orbital parameters we use to characterise a planet in- 
fluence is: the pulsar mass (Mpspr), the “projected mass” of the 
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Table 1. Table showing the period ranges and corresponding number of 
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samples. 
min Py [d] max Py [d] N 

21.3 42.5 2x 104 
42.5 85 1x 104 
85 170 8 x 103 
170 340 6x 10° 
340 390 6 x 103 
390 780 4x 103 
780 1560 4x10 
1560 3120 4x10 
3120 6240 4x10 


planet (m sini), the orbital period (Pẹ), the eccentricity (e), the ar- 
gument of periapsis (w) and the planet phase on its orbit at time 
55000 MJD (¢). 

We fix Mpsp at the representative value of 1.4 Mo (e.g. Lattimer 
2012), and sample over the other orbital and noise parameters using 
the EMcEE ensemble MCMC sampler (Foreman-Mackey et al. 2013, 
2019) within ENTERPRISE. Uniform priors are used for the orbital 
parameters w, e, ¢, allowing w between O and 27, e between 0 
and 0.9, @ between 0 and 1, and for each of the nine Py bins. For 
m sini, a log-prior is used as we want to explore several orders of 
magnitude between 1074 and 100M@. Note that w and @ can be 
highly correlated, as they both depend strongly on the position of the 
periapsis, which is poorly defined for low-eccentricity orbits. 

For each pulsar, the analysis is split into eight period ranges, uni- 
form in log-space, as well as a narrower range around Py = 1 year, 
since fitting for the pulsar position in the timing model can absorb 
power at the orbital period of the Earth. Table 1 shows the period 
ranges (min P,—max Py), as well as the corresponding initial number 
of samples per Markov chain used by the ensemble sampler (N). The 
number of samples for each period range was chosen as a compro- 
mise between the running time of the sampler and having a sufficient 
number of samples such that injected planet signals were consistently 
found. For each run of N samples, 64 chains are used, which after 
“burnin” and “thinning” results in 4.8N independent samples of the 
posterior. 

To check for 3-a detections and compute the 95-percentile mass 
limit, the mass posterior is re-weighted to reflect a uniform linear 
prior and detections are flagged if the mean is more than 3-0 away 
from zero. 

An exact treatment would require solving Eq. 4 for the value of 
E(t), at each ToA and each e. Due to the large amount of data to 
process, we instead interpolate over pre-computed values of E(t) over 
a fixed grid of 1,000 mean anomaly values and 90,000 eccentricity 
values, chosen such that the maximum fractional error is less than 
1078. This approximation reduces the run-time by a factor of ~ 6. 


5 RESULTS 
5.1 Relevant examples 


As the only pulsar known to host Earth-mass planets, PSR B1257+12 
provides a useful example for our analysis. Two of these planets are 
in the orbital period range of our search - one at Py = 66.5 days, and 
the other at Pp = 98.2 days (Konacki & Wolszczan 2003). Panel a 
of Fig. 1 shows the results of our pipeline for B1257+12 when no 
planetary companions were included in the initial pulsar timing pa- 
rameters. The two planets are easily detected at the expected orbital 
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Figure 1. Example plots of outputs for three different pulsars after our analy- 
sis. Panel a shows the results for B1257+12 without including its two known 
higher-mass planets of orbital periods 66.5 and 98.2 days in the timing model; 
panel b shows the results for JO134—2937, which gave the lowest mass limits 
for a normal pulsar; panel c shows the results for J1744—1134, which gave 
the lowest mass limits for a millisecond pulsar. In each plot, the step line 
shows the 95% mass limits, such that we consider masses above this line 
to be highly unlikely; the squares show the maximum-likelihood values of 
period and mass for each period bin. Furthermore, in panel a, the horizontal 
and vertical dotted lines show the known values of the corresponding orbital 
parameters for the two known planets, while the “plus” markers show 3-0 
detections as a result of our analysis. Note that the x-axis increases to the left. 


period and projected mass. Note that we have not included the up- 
per limits in the period bins in which a detection is made to avoid 
confusion. 

The constraints on the maximum mass of an undetected planet 
vary for each pulsar, depending on the length and cadence of the 
observations available, the presence of glitches or known binary 
companions, as well as on the amplitudes of any white and red 
noise present. The mass limits are generally better constrained for 
millisecond pulsars than for normal pulsars. This is illustrated in 
Fig. 1, which shows the lowest mass limits for a normal pulsar (PSR 
J0134—2937; panel b), and a millisecond pulsar (PSR J1744—1134; 
panel c), respectively. At the most extreme, we can put mass limits 
as low as 1/100th of the mass of the Moon (or ~ 10-4 Me). 


5.2 3-0 detections 


Excluding B1257+12, we found 15 pulsars in our sample for which 
the analysis indicates a greater than 3-o detection in at least one 
of the searched period bins. None of these are millisecond pulsars. 
An algorithmic detection on its own, however, is not necessarily evi- 
dence for an orbiting planet, and further inspection is warranted. The 
methodology described in Coles et al. (2011) was used to estimate 
the power spectra from our data for each of the pulsars showing de- 
tections. These were used to visually inspect the agreement between 
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Figure 2. The power spectral density (PSD) against frequency for PSR 
B1540—06 is shown by the black continuous line. The fitted power-law, with- 
out including a planet influence, and as given by Eq. 6, is shown by the dashed 
red line, with the value of the index y given in the legend. The white noise 
level is shown by the horizontal blue dashed line. The detected periodicity, as 
given in Table 2, is shown by the vertical green line, and a peak in the PSD 
can be seen at its approximate position. 


the noisy power spectra and the power-law model, as well as the be- 
haviour at the detected period. An example is given in Fig. 2, which 
shows the estimated power spectral density (PSD) as a function of 
frequency, as well as the fitted power-law model without including a 
planet influence. The equivalent plots for the other 14 pulsars can be 
seen in Appendix A. As many of these pulsars are dominated by red 
timing noise, we found it useful to additionally visualise residuals af- 
ter subtracting the lowest frequency components of the timing noise 
model. We limit this to frequencies below 1/2P, so that we avoid 
falsely emphasising the detected periodicity. Additionally, in order to 
better understand the significance of the planetary signal, for any bin 
in which such a detection was made, we repeated the extended EN- 
TERPRISE analysis using a nested sampler (pyNEsTy; Speagle 2020). 
This nested sampler was used to estimate the evidence and there- 
fore the Bayes factor between models with and without the planetary 
companion included. We consider a log-Bayes factor larger than ~ 3 
(i.e. odds ratio larger than ~ 20) to be “strong” evidence (Kass & 
Raftery 1995). We categorise and discuss these detections in the sec- 
tions below, as well as summarise the main parameters of the fitted 
planets found in our analysis in Table 2. 


5.2.1 Pulsars with known quasi-periodic spin-down 


Quasi-periodic (QP) timing irregularities have been observed in pul- 
sar data as variations in their frequency (v) or frequency derivative 
(v). Lyne et al. (2010) characterised several pulsars showing QP be- 
haviours, and concluded that these were caused by changes in the 
pulsars’ magnetosphere, as they were correlated with pulse profile 
changes. PSRs B1540—06, B1714—34 and B1826—17 were among 
those found by Lyne et al. (2010) to show QP timing noise, and 
the frequencies of their periodicities were estimated as 0.24(2) yr7!, 
0.26(4) yr~! and 0.33(2) yr~!, respectively, where the numbers in 
brackets represent the 1-o error in the last digit. 

Similarly, PSR B1828—11 is known to show strong QP oscillations 
in v, at a frequency of 0.73(2) yr~!. The second harmonic of this is 
also clearly detectable (Lyne et al. 2010; Stairs et al. 2019). In our 
analysis, we detected an oscillation consistent with a planet for each 
of PSRs B1540-06, B1714—34 and B1826-17, with orbital periods 
which correspond to frequencies indistinguishable from those mea- 
sured by Lyne et al. (2010). Note that the orbits found for these planets 


Table 2. Table showing the orbital period, mass and eccentricity values of the 
fitted planets found in our search. The values in brackets correspond to one 
standard deviation. The values of eccentricity given as an upper limit (e.g. < 
0.14) represent 95% upper limits of values consistent with 0. The last column 
shows which section in this paper discusses the respective pulsar. 


PSR name Py[d] msini[Me] eœ Sect. 
B1540—06 1473(14) 112) 0.12(4) 5.2.1 
B1714-34 1417(7) 6.3(6) 0.14(3) 5.2.1 
B1826-17 1102(5) 2.6(4) 0.35(6) 5.2.1 
B1828-11 (a) 231(6) 1.3(2) <0.14 5.2.1 
B1828-11 (b)  498(2) 6(1) 0.23(5) 5.2.1 
B0144+59 319(1) 0.060(8) <045 5.2.2 
B1727-33 350(1) 3.5(5) <0.26 5.2.3 
B2053+36 1013(9) 0.09(2) <040 5.2.3 
J1758-1931 719(8) 6.1(9) <043 5.2.3 
J1843-0744 650(7) 1.0(2) 0.4(1) 3.2.3 
J1904+0800 946(15) 1.0(2) < 0.18 be 
J2216+5759 117(9) 3.5(6) < 0.41 DAS 
J2007+3120 723(8) 2.3(3) < 0.38 5.2.4 
J1947+1957 1070(9) 3.7(3) 0.56(6) 5.2.5 
B1931+24 51801160) 56(8) 0.25(7) 5.2.6 
B0823+26 28.0(1) 0.08(2) 0.37(2) 5.2.7 


are eccentric, such that the oscillations in the residuals are asymmet- 
rical. For B1828—11, we find that, within our model, the residuals 
are best described by the influence of two planets, at the expected 
orbital periods corresponding to the first and second harmonics as 
estimated by Lyne et al. (2010). 

Fig. 3 shows the residuals for PSRs B1540-—06, B1714—34, 
B1826-—17 and B1828—11, respectively. Panel b in each pulsar figure, 
showing the residuals without the low-frequency red noise compo- 
nents, as well as the corresponding Rømer delay of the fitted planet(s), 
illustrates that the behaviour is not perfectly periodic, and the purely 
periodic orbital solution only removes a fraction of the power on these 
time scales. While a model fitting for one (two) periodic term(s) bet- 
ter describes our data than the model without it (as also supported by 
the log-Bayes factors, which are larger than 4 in all cases), a quasi- 
periodic term, caused by intrinsic properties of the pulsars, would be 
more appropriate, and would rule out the planetary-like detections 
for these pulsars. 


5.2.2 B0144+59 


PSR B0144+59 is dominated by red timing noise, but we detect 
a significant 319(1) day periodicity in our analysis. Removing the 
lowest frequency red-noise components reveals a highly periodic 
signal in the residuals, especially clear in the more recent DFB data 
where the precision of the observations is higher, as can be seen in 
Fig. 4. The existence of this periodic signal is also supported by the 
Bayesian evidence, as calculated using the software DyNESTY: our 
model has a log-Bayes factor of 9.2 + 0.5 compared to the standard 
ENTERPRISE model, which does not include a planet. 

Given the detections in pulsars with known correlations between 
profile shape and spin-down presented in Section 5.2.1, we looked 
for pulse shape changes in PSR B0144+59. As shown in Figure 5, 
PSR B0144+59 exhibits an interpulse, but only the main pulse is 
strong enough to detect observation-to-observation variations. We 
used PSRSALSA (Weltevrede 2016) to fit the components of the pulse 
profiles using von Mises functions and to estimate the full-width 
half-maximum (FWHM) of the main pulse at each observation. Fig. 6 
shows the measured width overlaid with the ý derived from our best- 
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Figure 3. The residuals of PSR B1540-—06, B1714—34, B1826-17 and B1828-—11 are shown. The top panel (a) shows the residuals after only removing the 
effect of deterministic timing parameters, as they are usually seen when using e.g. TEMPO2. The middle panel (b) shows the result of also removing the low 
frequency red-noise components, up to half the frequency of the fitted planet; the green solid line in the same panel is the Rømer delay that the fitted planet would 
cause. Finally, the bottom panel (c) is the result of removing the planet effect from the residuals in the middle panel. Note that the timing and noise parameters 
used here were those found by our analysis, which included a planetary companion. 
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Figure 4. The residuals of PSR B0144+59 are shown. For a detailed descrip- 
tion see the caption of Fig. 3. 


fit orbital model. The strong correlation, with correlation coefficient 
of ~ 0.7, is strong evidence that the observed periodic variations 
in PSR B0144+59 are magnetospheric in origin, similar to that of 
the pulsars in Lyne et al. (2010). Therefore, we conclude that the 
highly-periodic behaviour found in the residuals of PSR B0144+59 
is likely not due to an external orbiting body, but to an intrinsic 
magnetospheric effect of the pulsar. 

PSR B0144+59 exhibits very periodic spin-down oscillations, 
however this is not uncommon amongst the pulsars that exhibit QP be- 
haviours (see e.g. the pulsars in Section 5.2.1). Our analysis searches 
for purely periodic variability and therefore selects for those with 
the most stable periodicity. Nevertheless, this does emphasise the 
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Figure 5. The different main pulse profile shapes for BO144+59 are shown in 
the top panel, including an average profile in black. The average shape of the 
interpulse is also shown in the bottom panel, while there are large variations 
between different interpulse shapes due to the low signal-to-noise. 


somewhat surprising result that pulsars can exhibit extremely stable 
variability with periods of tens or hundreds of days. 


5.2.3 Additional quasi-periodic pulsars 


In this section, we discuss six of the pulsars for which we detected a 
planet-like periodicity, as they show similar behaviours, namely PSRs 
B1727-33, B2053+36, J1758—1931, J1843—0744, J1904+0800 and 
J2216+5759. 

Fig. 7 shows the residuals of these pulsars. The systematic be- 
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Figure 6. The estimated full-width-half-maximum (W50) values character- 
ising the pulse profile shape variation of PSR B0144+59 are shown by the 
black data points. The red line shows the ý estimated using our fitted planet 
model. 


haviours of the residuals after subtracting the planet-like influence 
(shown in panels c) suggest that these Rømer delays are not enough 
to account for the shape of the oscillating residuals of these pulsars. 
Furthermore, we do not expect these residuals to be caused by extra 
companions in the system, since their shape does not indicate a differ- 
ent periodic oscillation. Unlike PSR B0144+59, the signal-to-noise 
of the observations of these pulsars is insufficient for us to reliably 
search for pulse profile variations. 

PSR B1727-—33 has two known glitches, at observation times 
52 107 MJD and 55 930 MJD, respectively (Jodrell Bank Glitch Cat- 
alogue?; Espinoza et al. 2011). We include the standard glitch fit- 
ting parameters, as well as the TEMPO? glitch recovery parameters 
(GLFOD and GLTD), in the fit for our analysis. Note that, although 
the period of the fitted planet is close to one year (Pp = 350(1) d), we 
can rule out the effect of the rotation of the Earth as a source of this 
oscillation, as our analysis includes a fit for the position and proper 
motion of the pulsar simultaneously with all other parameters. 

We therefore find that the pulsars in this section seem more 
likely characterised by similar QP behaviours to the pulsars we dis- 
cussed in Section 5.2.1 than orbital motions, but without further 
data we cannot be certain of the origins of the detections in PSRs 
B1727-33, B2053+36, J1758—1931, J1843—0744, J1904+0800 and 
J2216+5759. 


5.2.4 J2007+3120 


Our initial analysis of PSR J2007+3120 revealed an oscillation con- 
sistent with a planetary companion of orbital period 723(8) days. The 
corresponding Rømer delay is shown in green in panel b of Fig. 8, 
on top of the residuals of this pulsar. After removing the influence of 
this planet, the residuals (shown in panel c of Fig. 8) further show a 
periodic behaviour. When repeating our analysis to fit for two planets 
accordingly, the parameters of the second potential companion were 
less constrained, giving an orbital period of 1297(76) days. 

The log-Bayes factors between models with zero, one and two 
planets suggest that, while the existence of the first planet (of shorter 
period) is strongly preferred (by 5.8 + 0.5), there is not a strong 
preference between a red noise component and the second planet, 
since the two models have indistinguishable evidences. 

While the residuals of J2007+3120 appear highly periodic, the 
time span of our current data only allows for a few oscillations. It 
follows that this is not sufficient to attest whether (any of) these two 
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planet-like influences are indeed due to planetary companions to the 
pulsar, or are due to other effects, such as intrinsic QP noise. The 
signal-to-noise in the observations of J2007+3120 is very low, which 
prevents us from reliably searching for pulse profile variations. 

The two periodicities found in our analysis are consistent, within 
two standard deviations, to a 2:1 ratio. On one hand, this ratio is char- 
acteristic of a harmonic relationship shown by QP noise, as is found 
for PSR B1828-—11. On the other hand, an orbital resonance with a 
2:1 period ratio has been found to be fairly common in extrasolar 
planetary systems (Lissauer et al. 2011). If the periodic behaviour is 
indeed due to planetary companions, detecting the effect of a (near) 
resonance in the form of a predictable variation of the orbital param- 
eters would firmly prove the existence of these companions around 
J2007+3120. This is the same technique used to confirm the first two 
planetary companions of B1257+12, which are also in resonance, 
with a period ratio of 3:2 (Rasio et al. 1992; Wolszczan 1994). Fol- 
lowing the three-body analysis of Malhotra (1993), we estimate that 
for the parameters of our two found planets around J2007+3120, the 
corresponding resonant gravitational interactions would introduce 
deviations in the residuals of the pulsar as an oscillation of period of 
a few tens of years and amplitude of ~ 1s. This is undetectable with 
our current data, as the residuals are characterised by a white noise 
level of ~ 1000us. However, long term observations using a highly 
sensitive instrument such as FAST, CHIME or the Square Kilometre 
Array might be able to detect/rule out these oscillations, and therefore 
settle the cause of the periodicities seen in J2007+3120. 


5.2.5 J1947+1957 


PSR J1947+1957 was discovered by Nice et al. (2013), who at the 
same time also claimed the detection of a glitch in the ToAs at 
MJD 55 085. In our analysis of J1947+1957, we found a planet-like 
influence, characterised by an orbital period of 3 years, projected 
mass of 4Mg@ and a highly eccentric orbit (e = 0.6); this model 
was preferred to the standard ENTERPRISE fitting by a large log-Bayes 
factor, of 7.9+0.5. 

The residuals of PSR J1947+1957 (see Fig. 9) are dominated by 
five peaked oscillations. The purely periodic model is successful 
at removing a good deal of the structure in the timing residuals, 
though the remaining residuals in panel c still have structures that 
appear correlated with each periastron passage. Although in theory 
the pulsar wind or tidal effects could perturb the orbit at periastron, it 
is worth noting that the closest approach of ~ 1 AU for the planetary 
companion is not especially small, and such effects could also be 
caused by variability in the periodicity of the oscillation. 

Alternatively, the cusps observed in the timing residuals could be 
due to glitches, as originally proposed by Nice et al. (2013) for the 
event at 55 085 MJD. There is no indication in our data of the presence 
of a glitch recovery for this proposed glitch. The subsequent four 
features are much less well sampled, so it is hard to tell if they exhibit 
the sudden change of gradient expected for a pulsar glitch. Although 
not common, some frequently glitching pulsars have been observed to 
exhibit glitches at regular intervals (e.g. Yu et al. 2013; Fuentes et al. 
2019; Basu et al. 2020). However, note that the inferred glitch size 
for J1947+1957 would be much smaller than those of other observed 
quasi-periodically glitching pulsars. Another potential origin of the 
periodic behaviour in J1947+1957 could be an intrinsic QP effect 
in the pulsar. The very low signal-to-noise of our observations of 
this pulsar makes any study of the pulse profile variability extremely 
difficult. 

It is notable that the argument of periapsis found for the planetary- 
like signal in pulsar data was w = (80 + 20) deg, while the high- 
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Figure 7. The residuals of PSR B1727—33, B2053+36, J1758—1931, J1843-—0744, J1904+0800 and J2216+5759 are shown. For a detailed description see the 


caption of Fig. 3. 


eccentricity, yet symmetrical shape of the cusps (as seen in Fig. 9) 
also suggests that the planetary orbit would have to have a special, 
symmetrical orientation in the sky, such that the observer, the peri- 
astron and the apastron would be aligned, i.e. w ~ 90 deg. However, 
this symmetry is typical in QP magnetospheric pulsars. 


5.2.6 B1931+24 


PSR B1931+24 was the first intermittent pulsar discovered; it is 
known to show normal pulsar behaviour for 5 — 10 days, and then 
undergo extreme nulling events for 25 — 35 days, making it unde- 
tectable (Kramer et al. 2006). The possibility that this observed QP 
behaviour is due to a binary companion of B1931+24 has been stud- 
ied by Rea et al. (2008) and Mottez et al. (2013); both concluded that 
the influence of a companion would not be enough to account for all 
the properties of the intermittent signal. 

In our analysis of this pulsar, we do not find that a single planet-like 
periodicity in the range 21.3 - 42.5 days is a good model to describe 
the residuals. This is because, as described by Young et al. (2013), 
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the nulling events do not have an exact periodicity, but multiple, 
narrowly-spaced periodicities in the range discussed. 


However, we find a different planet-like influence, of period ~ 
14 years. Fig. 10 shows the only two oscillations that correspond 
to our current data span. The shape of the residuals in panel b as 
compared to the Rgmer delay and the clearly systematic behaviour 
of the residuals in panel c of Fig. 10 suggest that this model does 
not fully describe the behaviour seen, and therefore this is likely not 
actually caused by a binary companion. Nonetheless, the log-Bayes 
factor between accounting for a planet-like periodicity, and not, is 
7.9 + 0.5. This relatively large value suggests that there should be a 
better model to describe the current residuals than the purely power- 
law description of the red noise used in this analysis. However, with 
only two oscillations seen in these residuals, we cannot determine if 
this pulsar will continue showing the same periodic behaviour, and 
therefore what the cause of this is. 
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Figure 8. The residuals of PSR J2007+3120 are shown. For a detailed de- 
scription see the caption of Fig. 3. A second planet-like periodicity fit is 
shown by the solid line in panel c. 
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Figure 9. The residuals of PSR J1947+1957 are shown. The dotted vertical 
line corresponds to the glitch found by Nice et al. (2013). For a detailed 
description see the caption of Fig. 3. 
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Figure 10. The residuals of PSR B1931+24 are shown. For a detailed de- 
scription see the caption of Fig. 3. 


5.2.7 BO823+26 


PSR B0823+26 is one of the brightest known radio pulsars in the 
northern sky, and has been observed for nearly 50 years. This pulsar 
is known to exhibit emission phenomena such as subpulse drifting 
and mode-changing (Sobey et al. 2015; Basu & Mitra 2019), as well 
as nulling over a range of timescales from minutes to hours, often in 
clustered groups and with a fraction of ~ 7% (Herfindal & Rankin 
2009; Redman & Rankin 2009; Young et al. 2012). 

Our initial analysis flagged a planet-like signature in the residuals 
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Figure 11. Summary of 95-percentile mass limits for all pulsars in the anal- 
ysis. Note that the x-axis increases to the left. The colours used represent 
the fraction of pulsar mass limits lower than the respective mass, at each 
corresponding period bin. The diamonds represent the two known planets of 
B1257+12 whose orbits are within our parameter space, for reference. The 
stars show the detections flagged by this analysis, excluding the known QP 
pulsars presented in Lyne et al. (2010). 


of B0823+26, with a period of 28(2) days. The residuals of this 
pulsar show that the influence of the fitted planet would be very 
small compared to the general white noise in the pulsar timing. On 
further inspection, the power spectrum (shown in Fig. A1) of the 
timing data for this pulsar appears more like a broken power-law 
with a steep red component transitioning to a flatter red component 
at a frequency close to that found by our analysis; this is likely to be 
a consequence of the complicated nature of its emission. The power 
spectrum of the noise from B0823+26 was therefore fitted poorly by 
a single power-law model, resulting in a false planet detection. 


5.3 Discussion on the results for the full data set 


We show a summary of the 95% mass limits of all the pulsars in 
our sample in Fig. 11. The different contours illustrate the fraction 
of the pulsars in our sample that allow a planet below the respective 
projected mass within the 95% limit, for each period bin. Note that 
we can only estimate the projected mass on the sky with respect to 
the observer (m sini); therefore the actual physical mass limits are 
on average ~ 20% larger, assuming a random orientation of binaries 
in the sky (since the mean of sin i is 0.785). 

In general, the mass limits as presented in Fig. 11 depend on the 
properties of the pulsars, as well as on the properties of the obser- 
vations. Mass limits close to a 1-year orbital period are noticeably 
worse than all the rest due to the required fitting for the proper motion 
and distance to the pulsar, which vary periodically with the Earth’s 
orbit and therefore reduce our sensitivity. Similarly, any known bi- 
nary companions reduce the sensitivity of our search in the respective 
period bin. At the long-period end the sensitivity is mainly limited 
by the amount of red noise in the residuals. At periods larger than 
5-10 years, the total timespan of the observations can also limit the 
sensitivity of our search, especially for more recently discovered pul- 
sars. In the short-period limit the white noise is the dominant factor 
in limiting our sensitivity. 

Fig. 11 also shows a summary of the planet detections found in 
our analysis, with the exception of the known QP pulsars included 
in Lyne et al. (2010). As discussed in the preceding sections, we 
suspect that the majority of these are likely to be spurious detections 
caused by QP noise in the pulsar residuals. Interestingly, these show 
up in a pure periodicity search, emphasising that pulsars can exhibit 
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Figure 12. Joint period-mass posterior for all pulsars. The red, green, orange 
and blue lines show the 95%, 68%, 32% and 5% mass limits, respectively, 
of the mass posterior formed by adding the mass posteriors of all the pul- 
sars, except those that showed a detection. This is done separately for each 
period bin, as can be seen. The number shown in the caption (e.g. “i.e. 740”) 
represents the number of pulsars corresponding to the respective percentage 
(e.g. 95%). Two curves from K15 are also included for comparison, showing 
the 95% mass limits for the pulsars in their sample of characteristic ages 
above and below 0.3 Myr. The grey dots represent known exoplanets, from 
the Exoplanet Database (Wright et al. 2011). Note that the x-axis increases to 
the left. 


highly stable periodic variability on timescales of hundreds of days. 
A detailed re-analysis of the sources in Lyne et al. (2010) is ongoing 
(Shaw et al., in prep). 

It is notable that in all cases where the eccentricity is well con- 
strained the eccentricities are between 0.1 and 0.6, and are certainly 
not the circular orbits that might be expected from formation in a 
fall-back debris disk. This further supports the hypothesis that these 
detections are more likely due to intrinsic spin variation rather than 
orbital motion. 

In the period range of 0.2 — 4.3 years, with the exception of the 
aforementioned small interval around 1 year, about 70% of the pulsars 
in our sample have mass limits lower than 1 Mg. This implies that 
the majority of our pulsars are generally unlikely to host any planets 
larger than this. Furthermore, 5-10% of all the pulsars in our data set 
have mass limits lower than one Moon mass (or ~ 0.01 Mẹ ) in the 
same period range, ruling out any substantial planetary companions 
in our search range. Note that these mass limits are similar to the total 
mass of the asteroid belt model around PSR B1937+21 as estimated 
by Shannon et al. (2013). The largest known mass of an asteroid 
is, however, two orders of magnitude smaller (~ 1074 Mọ; Baer & 
Chesley 2008). 

Fig. 12 shows another way of illustrating the mass limits for our 
entire pulsar data set. This is an approach similar to that of Fig. 4 of 
K15, and we include some of their results for comparison. To obtain 
the curves shown in Fig. 12, the mass posteriors of all the pulsars in 
our data set, except those that showed detections, were added for each 
period bin to create a mass posterior of the data set. The 95%, 68%, 
32% and 5% limits were then calculated and are shown by the red, 
green, orange and blue lines in Fig. 12. Note that our mass sampling 
is re-weighted into a uniform prior between 1074 and 100Mg. 

Comparing our results with the limits from K15, we note that 
our data set includes more than 5 times the number of pulsars in 
the K15 sample (consisting of 151 pulsars). Therefore, although for 
most period bins the 95% mass limits of K15 give better constrains 
than our 95% mass limits, they represent a much smaller number 
of pulsars, i.e. 50-100 for K15 as opposed to 740 for this analysis. 
We also note that K15 specifically used pulsars of high spin-down 
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luminosity and with regular cadence observations, whereas our sam- 
ple consists of essentially all pulsars observed by the Jodrell Bank 
timing programme, and includes a number of pulsars which have 
been infrequently and irregularly observed. 

The mass limits shown can be interpreted as the fraction of pulsars 
in the data set that do not show a detection below the respective 
projected mass. For example, we can conclude that roughly two- 
thirds (or 68%) of the pulsars in our data set (i.e. ~ 530 pulsars) 
do not show any planets above 2 Mg and with orbital period smaller 
than 1 year, as well as any planets above 8 Mg and with periods larger 
than 1 year. Likewise, for about 40 of our pulsars, we can rule out 
planetary companions of projected masses above 0.2 Mg. Despite 
the large increase in sampled pulsars, our mass limits are better than 
those of K15 for large orbital periods. This is a consequence of 
our longer-span data sets, which improve our potential for detecting 
longer-period orbits, and often allow us to better distinguish between 
red noise processes and truly periodic behaviours. 

This large sample study shows that although planets around pulsars 
are rare, there is still a large population of pulsars for which we cannot 
rule out planets of order ~ 1 Mg. The prevalence of QP timing noise 
in pulsars exacerbates this problem as it can lead to spurious signals 
that greatly increase the detection threshold. Nevertheless, we can 
certainly rule out a population of planets greater than ~ 10 Mẹ from 
our pulsars, even those with extensive timing noise. This also rules 
out a population of Jupiter-mass companions (including so-called 
“diamond planets”) with orbital periods in our search range. 

As seen in Fig. 12, most of the known exoplanets could be ruled 
out around our sample of pulsars. Recent results suggest that the 
distribution of exoplanet-to-host-star mass ratio (q) follows a broken 
power-law, with maxima at q ~ 3 x 1075 for systems with host stars 
being “G dwarfs” (Pascucci et al. 2018), and at q ~ 2 x 1074 for 
“microlensing” systems (Suzuki et al. 2016; Udalski et al. 2018), 
respectively. One of the pulsar planet formation models requires the 
planet(s) to have initially formed around a massive star and survive 
the subsequent supernova state (see e.g. Bailes et al. 1991; Podsiad- 
lowski 1993). If we were to assume the shape of the microlensing 
q-distribution holds when extrapolated to host stars massive enough 
to be pulsar progenitors, i.e. more than ~ 8 Mg (see e.g. Woosley 
& Weaver 1986), the most likely planets would have masses much 
above 100 Mg, and anything below 10 Mg would be extremely un- 
likely. Therefore, as also argued by Martin et al. (2016), our analysis 
would support the idea that fossil (or remnant) planets are, at most, 
rare. If, on the other hand, we were to make the even more speculative 
assumption that planets forming around pulsars would follow similar 
q-distributions as more typical star systems, we would also expect 
that these exoplanets were 10 — 100 Mg, which are absent from our 
search. 

Alternatively we might consider our results in the context of PSR 
B1257+12, which is the only pulsar known to host terrestrial (i.e. 
~ Earth-mass) planets. The joint posterior shown in Fig. 12 implies 
that the probability that a system with a planet of mass similar to 
that of the larger planets of B1257+12 (~ 4Mg@) would have gone 
undetected in our sample of pulsars is less than 0.3. This implies 
a reasonable upper bound of 3 undetected B1257+12-like planetary 
companions amongst our 800 pulsars, with a probability of less than 
0.03. Given our broad sample, this implies that no more than 0.5% 
of pulsars host terrestrial planets as large as those in PSR B1257+12. 
We therefore confirm the hypothesis that the formation of planets 
around pulsars is rare, and PSR B1257+12 is a special case. It is 
worth noting that the smaller planet in PSR B1257+12 (0.02 Mẹ) 
would be undetectable in 95% of our sample, and we cannot rule 
out a large population of tiny planets or asteroids. However, given 
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the unique case of PSR B1257+12, it is unclear if such tiny planets 
are likely to form without the presence of larger planets. Further, 
such planets may be part of larger debris disks or asteroid belts and 
exhibit noise-like timing behaviour, rather than the purely periodic 
signal considered in this work (Jennings et al. 2020). Further study 
into the existence of a population of these systems with only small 
planets would help determine how common a system such as that of 
B1257+12 truly is. 


6 CONCLUSIONS 


In this work, we have presented the largest search for planets or- 
biting pulsars to date, analysing JBO observations of 800 pulsars. 
We confirm that PSR B1257+12 must have an unusual formation 
mechanism, placing an upper bound of 0.5% of pulsars exhibiting 
similar planets. We rule out a population of undetected planetary 
companions greater than ~ 10 Mg. The timing noise present in most 
pulsars means that we cannot rule out a substantial population of tiny 
(< 0.1 Mg) planets, though it is not clear if such planets would exist 
in isolation. 

For most pulsars in our data set, no planet-like periodicity was de- 
tected. Of the few that did show these detections, several are known 
to exhibit QP spin-noise which is correlated with variations in their 
pulse profile. We have identified that PSR BO144+59 also exhibits 
this behaviour, by looking at the variation of its pulse profile shape. 
We suspect that many of the detections made in our analysis are 
caused by similar effects, but the signal-to-noise is insufficient to see 
pulse profile changes. The prevalence of QP noise in pulsar resid- 
uals makes searches for low-mass planetary companions to pulsars 
very challenging and further work is needed for robust pulsar timing 
analysis in the presence of such noise. 
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APPENDIX A: POWER SPECTRAL DENSITY PLOTS 


The plots of power spectral density against frequency for the pulsars 
showing detections are shown in Fig. Al (except PSR B1540—06, see 
Fig. 2). The frequency range shown is limited at the low end by the 
timespan of our observations, and at the high end by our chosen value 
of 0.1d7!. Notably, the intermittent PSRs B1931+24 and B0823+26 
show that the model described by a single red noise power-law and a 
flat white noise component does not describe the data well. Specifi- 
cally in the case of B0823+26, the extra, unmodelled power at high 
frequencies leads to what we consider a false periodicity detection at 
28 days. 

The remaining plots show that, in general, our periodicity detec- 
tions correspond to peaks in the power spectra, unmodelled by the 
red noise power-law, as expected. However, note that this power-law 
model is not a perfect description of the variations in the power spec- 
tral density, and therefore a more involved analysis is necessary, as 
described in the main body of this work. 
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Figure A1. The power spectral density against frequency for the pulsars that showed detections, as summarised in Table 2. For more details see the caption of 


Fig. 2 
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